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We present an implementation of the Galerkin-Collocation method to determine the initial data 
for non-rotating distorted three dimensional black holes in the inversion and puncture schemes. 

The numerical method combines the key features of the Galerkin and Collocation methods which 
produces accurate initial data. We evaluated the ADM mass of the initial data sets, and we have 
provided the angular structure of the gravitational wave distribution at the initial hypersurface by 
evaluating the scalar ^'4 for asymptotic observers. 


I. INTRODUCTION 

The full 3D evolution of the Einstein field equations 
figures as the most-challenging task for numerical rela¬ 
tivity despite the progress achieved so far [l|. In order 
to evolve any 3D code, one needs to specify the initial 
data representing a physically relevant system. Among 
all possible configurations, those involving black holes are 
of interest. The strong gravitational fields produce the 
ideal arena in which the fully general relativistic effects 
take place. 

In this direction, there is a class of the black hole initial 
data known as distorted black holes that Bernstein et al 
0 introduced. They assumed the axisymmetric initial 
that consists in a black hole with or without rotation in 
interaction with a cloud of gravitational waves of variable 
intensity about the black hole. Later, Brandt et al @ 
relaxed the axisymmetry and considered the most general 
three-dimensional distorted black holes. An important 
motivation in establishing and evolving distorted black 
holes is to reproduce the late stages of binary black hole 
coalescence. In addition, the dynamics of distorted black 
holes can provide a simple framework to study in details 
the efficiency of gravitational wave extraction, together 
with the determination of wave templates perceived by a 
distant observer. 

Recently, we have applied the Galerkin-Collocation 
spectral method [J to determine accurately two initial 
data sets for numerical relativity: pure Brill waves, and 
axisymmetric non-rotating distorted black holes. These 
problems were considered previously in the realm of tra¬ 
ditional pseudospectral @ and finite difference methods 
0. Several relevant works dealing with pseudospectral 
codes for the determination of single black hole initial 
data can be found in Refs. M- 
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There are two main strategies to describe the initial 
data sets for single and multiple black holes. We mention 
the use of isometry conditions at the inner boundaries in 
order to represent the black holes throats 0. Another 
approach is the puncture method proposed by Brandt 
and Brugmann [l2|. This method proved to be very ef¬ 
fective in describing initial data for multiple black holes, 
in particular binary black hole systems. Basically, it con¬ 
sists in splitting the conformal factor of the spatial metric 
into singular and non-singular terms. Brown and Lowe 
[isjl applied the puncture method for the determination 
of distorted black holes spacetimes with the implemen¬ 
tation of adaptive mesh refinement in their finite differ¬ 
ence code to solve the elliptic equation resulting from 
the Hamiltonian constraint. As they have pointed out, it 
was necessary to perform the computation on a large grid 
with high resolution near the black hole. The puncture 
data for binary black holes or neutron stars using single 
and multi-domain s pec tral methods were considered by 
Grandclement et al ll J. Ansorg et al fsL [TR , Pfeiffer 
[IiEIj Foucart et al |l8j|. Ruchlin at al [T^, Lovelace et 
al [23 and Koutarou at al [2l|. 

The main goal of the present work is to apply the 
Galerkin-Collocation method [1, [2^ to obtain three- 

dimensional distorted black holes initial data sets. We 
have developed algorithms in the realm of the inversion 
method and the puncture method with domain decom¬ 
position. We have organized the paper as follows. Sec¬ 
tion II presents briefly the basic equations of the 3 -I- I 
formulation for the initial data problem in both inver¬ 
sion and puncture methods. Section III is devoted the 
describe the numerical implementation of the Galerkin- 
Collocation method in both methods, where the choice 
of basis functions that satisfy the boundary conditions 
constitutes the cornestone of the codes. The spherical 
harmonics are the most natural basis functions for the 
angular domain, whereas the radial basis functions are 
expressed as suitable linear combinations of the Cheby- 
shev polynomials. The condition of inversion symmetry 
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had to be satisfied by imposing a relation between the 
unknown modes. We have implemented the puncture 
method with a simple version of the domain decompo¬ 
sition that consists in dividing the spatial domain into 
several regions, wherein each region we solve the Hamilto¬ 
nian constraint and match these solutions across the do¬ 
mains. Section IV shows the convergence tests together 
with the asymptotic behavior of the spin-weighted scalar 
which provides the pattern of the gravitational field 
associated to the distorted black hole at the initial slice. 
Finally, in Section V we make some concluding remarks. 


The metric of the initial hypersurface corresponding to 
a three-dimensional distorted black hole is expressed as 

ii, 

ds^ = ^-4 ^2 gjj^2 g,^^2j ^ ^5^ 

where the function q = q{r, 0, 4>) represents the distri¬ 
bution of gravitational wave amplitude that satisfy 
certain boundary conditions to ensure the regularity and 
the asymptotic flatness of the metric. These boundary 
conditions are: 


II. THE INITIAL DATA PROBLEM: BASIC 
EQUATIONS 


g(r,0,(/)) = q(r, 7 r/ 2 ,(/)) = 0, \mi q = 0{r ( 6 ) 

r^oo 


The basic equations for the initial data we are going 
to solve arise from the 3-1-1 formulation [ 2 ^ of the 
Einstein’s field equations. The initial data cannot be 
specified arbitrarily, but it must satisfy in vacuum four 
constraint equations given by. 


A + K‘^ - KijK^^ = 0 

(1) 

= 0, 

(2) 


where and Kij are the metric and extrinsic curvature 
of the 3-dimensional spacelike hypersurfaces that foliate 
the spacetime, respectively. All quantities are evaluated 
on the 3-dimensional hypersurfaces, and K = 

These four equations are known as the Hamiltonian and 
momentum constraints, respectively. 

Following Brandt et al Q we are going to consider the 
initial data at the moment of time symmetry, meaning 
that the extrinsic curvature is zero, Kij = 0 at the ini¬ 
tial slice or hypersurface. In this case, three constraint 
equations vanish identically remaining the Hamiltonian 
constraint — 0 which fixes the three-metric or initial 
data. It is appropriate to follow the York-Lichnerowicz 
approach expressing the metric 7 ^ in conformal 
form, 


7y — ^ Tijj (3) 

where 'k is the conformal factor and the metric jij are 
given. The Hamiltonian constraint becomes, 

--m = 0. (4) 

Here, and R are the Laplace operator and the Ricci 
scalar associated to the metric jij, respectively. There¬ 
fore, the Hamiltonian constraint © becomes an elliptic 
equation for the conformal factor 'k whose solution de¬ 
termines the initial data or the initial metric 7 ^j. 


We have considered the gravitational wave amplitude dis¬ 
tribution function introduced by to Bernstein et al 0, 


g(r, 0, (j)) = Aq sin" 6 




(I -|- CCOS^ (j)), 

(7) 


where Aq denotes the amplitude of the Brill wave [^ . 
the free parameter c indicates the deviation from axisym- 
metry, rj — ln(r/a) and n > 2 is an even integer; 770 , cr are 
constants associated to the position and width, respec¬ 
tively, of the Brill wave. 

The conformal factor must satisfy the condition, 


vk = l + ^ + 0(r-2) (8) 

at r —>■ 00 as a consequence of the Robin boundary condi¬ 
tion asymptotically, and parameter A4 is the ADM mass. 
The inner boundary is placed at the throat of the black 
hole and satisfies the isometry condition condition. 


/ 4' \ 

V^r 2a) 


= 0 , 


(9) 


where a = Mo/ 2 , and Mq is the mass of the black hole 
that results from setting q = 0. Therefore, the Bernstein 
data sets are obtained after solving the Hamiltonian con¬ 
straint in the region r > a satisfying the boundary con¬ 
ditions dH) and dH). 

An alternative way of constructing isometric distorted 
black hole data sets is provided by the so-called puncture 
method [H. The central idea of the puncture method 
is to split the conformal factor into a singular and non¬ 
singular terms, and to consider the whole spatial domain 
instead of being restricted to the region outside the throat 
of the hole. Accordingly, the conformal factor is written 
as, 


^' = u -I- 


m 


( 10 ) 
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where m is a new parameter of the method. Notice that 
the term TO/2r is singular at the origin, whereas the func¬ 
tion u is the nonsingular term. 

We present the Hamiltonian constraint expressed in 
function of u after substituting the decomposition m 
into the Eq. o, which results in, 


The spherical harmonics are complex functions imply¬ 
ing that the modes Ckim must be complex, but they sat¬ 
isfy certain symmetry relations in order to produce a real 
conformal factor. These symmetry relations are, 

Ch-m = (15) 


y^u--Ru=—R. (11) 

8 16 ^ ^ 

This equation is identical to Eq. with an addition 
source term \Qm/R. Brown and Lowe [l3j | have shown 
that the data sets obtained after solving Eq. m satisfy 
automatically the isometry condition if m = Mq = 2a. 
Eurthermore, it is possible to generate initial spacetimes 
that do not satisfy the isometry condition by setting m ^ 
Mo. 


since = (—1) ^Yim- As a consequence, the total 
number of independent unknown coefficients is (Nx + l) x 

(IV.+ 1)^- 

The inversion symmetry condition (IH]) is satisfied in an 
approximate way according to. 



(16) 


III. THE NUMERICAL SCHEME 
A. The inversion method 


for all I = 0,1, ..,1V, m = The integrals are 

evaluated using quadrature formulas, which is typical of 
the G-NI (Galerkin with numerical integration) method 
0: 


We follow the procedure we have employed in Ref. Q 
to deal with the axisymmetric distorted black hole data 
sets. The starting point is to establish an approximate 
expression for the conformal factor given by. 




j,k=0 

YL.{dj,<j>k)wjVk = 0, 


(17) 


N^,Ny I 

'i’a{r,9,4>) = 1+ ^ ^ CUm Xk{r)Yim{6,(t>), (12) 

k,l—0 m— — l 


where {9j,(j)k) with j = 0,1,.., iVi, fc = 0,1,.., N 2 , respec¬ 
tively, are the collocation points on the angular domain 
and given by. 


where Ckim represents the unknown modes and , Ny 
are the truncation orders that limit the number of terms 
in the above expansion. The angular patch has the spher¬ 
ical harmonics, Yim{9, (^), as the basis functions, whereas 
the radial basis functions, Xfe(r), are the same used for 
the axisymmetric case, 

Xk{r) = i(TLfc+i(r) - TLk{r)), (13) 

where, the rational Chebyshev polynomials TLk{r) are 
defined according to 

TLkir) =Tk(x =-— ^ (14) 

\ r — a + Lr J 

where Tk{x) is the traditional Ghebyshev polynomials. 
The radial domain a < r < 00 is equivalent to — 1 < 
X < 1 (see Eig. 1) with being the map parameter 
whose convenient choice improves the accuracy of the 
approximate solution. It can be shown that Xfc(r) = 
as r —^ 00 reproduces the boundary condition 
'I'(r, 0) = 1 + 0{r~^). 


Ink 

^ N 2 + I' 


The quantities Wj , Vk are the corresponding weights , 
and we have chosen Ni = N 2 = 2Ny for better accuracy 
in the calculation of the integrals (see also Ref. [2^.[30l|l. 
We have obtained {Ny + 1)^ linear algebraic relations 
for the unknown coefficients Ckim from (Ull) that can be 
solved to express the {Ny -|-1)^ coefficients CN^im in func¬ 
tion of the remaining ones. Therefore, we ended up with 
a total of Nx x {Ny + 1)^ independents coefficients. The 
use of the new coordinates x, y together with (j) covers 
the entire spatial domain as illustrated in Fig. 1. 

The residual equation associated to the Hamiltonian 
constraint is obtained by substituting the approximate 
conformal factor 'i’{r,9,(j)), into Eq. ([4]). We represent 
this equation by Res(r, 0,())) recognizing that it does not 
vanish identically due to the approximated conformal fac¬ 
tor. Next, the coefficients are determined such as to force 
the residual equation to zero in an approximate sense 311 
as we have done with the inversion symmetry condition. 
Then, it follows that. 
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FIG. 1. The angular variables {6, (j>) cover the surface of the sphere, or the subdomain r = constant (figure on the top). We have 
covered the entire spatial domain with the coordinates {x,y,(j>). We have shown in the left and right panels the subdomains 
X = constant (r = constant) and (j) = constant, respectively. Notice that the spatial infinity (r —^ oo) is placed at a; = 1. 
The initial data has equatorial plane symmetry due to the form of the gravitational wave amplitude with n even (see Eq. (0) 
(0 < 0 < n/2 or —1 < y < !)■ In addition, the angular dependence on rj) allows us to consider only the patch 0 < cp < <^/2. 


(Res, ^k{r)Yi^) = ^k{r)dr X (20) 

J a 

[ Res{r,e,(j))Yi*^{e,(j))dn = {), 

Jo. 

' -V-' 

Resim(j") 

where the (r) are known as the test functions [3l| , and 
I = 0,1, Ny, m = —l,—l + l, ..,1 and fc = I, 2,.., N^- As 
before we have evaluated the integrals over the angular 
domain using the Gauss quadrature formulas, or. 


of unknown coefficients Ckim- The preconditioning tech¬ 
nique (see Ref. Q and references therein) allows to re¬ 
duce the number of iterations in solving ill-conditioned 
linear systems especially those with an enormous num¬ 
ber of equations. In the present case we have reduced 
further the number of independent coefficients by taking 
into account the symmetries of the gravitational wave 
amplitude (cf. Eq. ([7])). Therefore, the resulting equa¬ 
tions are solved using standard linear solvers of Maple 
or Matlab, determining the coefficients and consequently 
the approximate conformal factor. 


Ni,N2 

Resimir)^ ^ Res{r,6j,(l)k)Y*^{9j,(j)k)wjVk. ( 21 ) 

j,k=0 

Now, by choosing test functions as delta of Dirac func¬ 
tions ^fe(r) = S{r — Tk), where rk represents the collo¬ 
cation points on the radial patch, we obtain the set of 
equations for the independent coefficients Ckim expressed 
as, 


B. The puncture method 

For the Galerkin-Collocation implementation of punc¬ 
ture method it is likewise necessary to establish an ap¬ 
proximate expression for the function u(r, 9, </>). The ful¬ 
fillment of the Robin condition implies that asymptoti¬ 
cally u = 1 + 0{r~^), and due to the nonsingular nature 
of M, it follows that, 


(Res, ^k{r)Yijn) = Resimirk) = 0, (22) 

where k = The radial collocation points 

are, 


— CL Lf 


(1 +Xk) 

1-Xk ' 



where 


(23) 


JVx.Afs I 

Ua{r,9,cj)) = l + E E Ckim Xk{r)Yi^{9,(j)). (24) 

k,l—0 m— — l 

This expression is identical to the approximate conformal 
factor given by Eq. (HI, where again the Ckim represents 
the unknown modes, Nx,Ny are the truncation orders, 
and the spherical harmonics are the angular basis func¬ 
tions. The radial basis functions are given by Eq. m, 
but the rational Chebysehv polynomials are given by, 


The set of equations (1771) has x {Ny -|- 1)^ linear and 
ill-conditioned algebraic equations for the same number 


TLk{r) 


Tk 



r-Lr \ 
r + Lr ) 


(25) 
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FIG. 2. Illustration of the regions X>i and II 2 viewed from 
the plane {r,9), where the bold line represents the boundary 
r = ro. 


in order to cover the whole radial domain 0 < r < 00 
being equivalent to — 1 < a; < 1. 

The determination of the coefficients Ckim follows the 
same steps we have devised previously. Since we are in¬ 
terested on the isometric data sets, we have set m = 2a 

A- . nr^ 

The domain decomposition technique [9|,l27| can be im¬ 
plemented more naturally in the scheme of the puncture 
method. It consists in dividing the spatial domain into 
two or more distinct regions, each one with approximate 
expressions for the function u. We present here a simple, 
but efficient version of the domain decomposition by es¬ 
tablishing two regions: the first region : 0 < r < rg 
and the second region defined by 2?2 : > ?'o- In Fig. 

2 we present the general scheme of the domain decom¬ 
position viewed in the plane (r, 0), where the bold line 
corresponds to the boundary r = rg separating the re¬ 
gions Vi and 'D 2 . Naturally, the junction conditions, 






dr 


ro 


du'd) 

dr 


(26) 


must be satified at the boundary r = r^. Numerically 
these relations are approximated, in the same way, as 
described by Eqs. (HU) and dni)- In Table 1, we summa¬ 
rize the approximate functions u at these regions in which 
the radial basis functions are chosen conveniently in each 
region, whereas the angular basis functions remain the 
same. These expansions have the spherical harmonics 
as the angular basis functions with the same truncation 
order Ny, but the truncation orders of the radial sector 
may be different in each region. 


IV. NUMERICAL RESULTS: CONVERGENCE 
TESTS, ADM MASS AND THE ANGULAR 
PATTERN OF GRAVITATIONAL RADIATION 

We present compelling numerical experiments showing 
the exponential convergence of the Galerkin-Collocation 
implementation for solving the initial data problem of 
three-dimensional distorted black holes using inversion 
and puncture methods. We have selected three tests in 
this direction: the convergence of the ADM mass, the 
L2-norm associated to the difference of the solutions cor¬ 
responding to successive truncation orders, and the L 2 - 
norm associated to the residual Hamiltonian constraint 
equation considering increasing radial resolution. 

We start with the calculation of the ADM masses of 
distorted black holes which are evaluated more efficiently 
using a formula derived by 6 Murchada and York [s^ . 

E-E = -^j^ (27) 

where A is the total energy of the hypersurface while 
is the energy associated to the conformal metric. As 
pointed out by Bernstein et al 0 this last term van¬ 
ishes since the conformal factor decays more rapidly 
than 1/r, therefore, the ADM mass is given by the in¬ 
tegral on the rhs of the above equation. By inserting 
Vq.'F = (d'l^/dr,l/rd'l^/d0,l/(rsm0)d'l'/dd) the final 
expression for the ADM mass becomes, 

Mavm = - ^lii^ ^ I dydd. (28) 

We have evaluated the above limit without approximat¬ 
ing the infinity to some finite radius r = Tmax- This 
feature is a consequence of defining the conformal factor 
in the whole spatial domain. Therefore, after obtaining 
the approximate conformal factor the ADM mass could 
be calculated by direct integration. 

We have obtained the convergence of the ADM mass by 
calculating the difference of the ADM masses correspond¬ 
ing to approximate solutions with distinct truncation or¬ 
ders. To be more specific, we have fixed Ny = 6 and es¬ 
tablished that 5M(Nx) = Madm(Nx + 7>) — Madm(Nx)- 
In Fig. 3 we present the exponential decay of 5M for 
the inversion method and the puncture method with do¬ 
main decomposition in which we have fixed = 30 
in the hrst domain. In both cases the saturation oc¬ 
curs at approximately = 70 (cf. Fig. 3), and 
Ag = cr = 770 = c= l,n = 4, but the puncture method 
with domain decomposition presents a better conver¬ 
gence rate. In the case of the puncture method, only 
the conformal factor defined in the second region, r > a, 
is used to calculate the ADM mass. The values of the 
ADM masses evaluated in both methods is unaffected 
to almost all significant digits. As a last comment, we 
have found that the puncture method without domain 
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Di : 0 < r < ro 


D 2 ■■ r > To 


ui^\r,e,cl,) ^ 1 + J2 ciZTLk{r)Yi^{e,^) 
Radial basis function: 

TL,{r) ^n{x=^-l) 


Ua^ {r,d,(t)) = 1 + J2 ‘^kfmXk (r}Yi^ ( 0 , (j>) 

Radial basis function: 

Xk{r) = ^{TLk+i{r) - TLk{r)) 
TLkjr) ^n{x= 


TABLE I. Approximate functions u(r, 6, cj)) defined at the regions Di and T> 2 - Notice that the radial basis functions are distinct 
in each region according to the boundary conditions. In the first region, we have considered the redefined Chebyshev polynomial 
functions with the map x = 2r/ro — 1 that connects 0<r<roto—l<a:<l. In the second region, we have used the same 
radial basis function defined for the inversion method. Although not indicated explicitly, we have adopted and as 
the radial truncation orders for the radial expansion in the first and second regions, respectively. 



FIG. 3. Exponential decay of SM{Nx) = Madm{Nx + 5) — 
Madm(Nx) for the inversion method (boxes) and puncture 
method with domain decomposition (circles) with = 30. 
In both cases Ny = 6 is fixed, as well Aq = 1, cr = 770 = c = 1. 
We have set Lr = 9.0 in both methods. 


decomposition in the realm of the Galerkin-Collocation 
implementation is not efficient in the sense of producing 
a poor convergence rate. 

We have noticed that in this first round of numerical 
experiments, the rate of convergence of the ADM mass 
is sensitive to the choice of the map parameter Ly. In 
both inversion and puncture method, we have found that 
Ly = 9.0 is the best value and adopted hereafter. The 
main criterion for choosing the map parameter is to co¬ 
incide it approximately with the scale of the problem 
under consideration (cf. Ref. pg. 369), but some 
trial-and-error was inevitable. 

For the convergence of the L2-norm of L2(<5'I'), we 
have considered the approximate solutions as described 
above. The calculation of L2(^^) for both inversion and 
puncture method takes into account the region outside 
the throat r > a. This quantity is given by. 


L2{6'S) 


I 

47r 



dy 


1/2 


(29) 

where 54'(a:, y, (j)) is obtained from ([8|) after changing the 
variables {r,9) to {x,y) according to 9 = arccosy and 
r = a + Lr{\ + x)/{I — x). We have used the property 
'I'(a;, y, (/>) = '^{x,—y,(j)) as a consequence of the reflec¬ 
tion symmetry about the plane 9 = 7r/2 or y = 0. The L2 
norm was calculated using quadrature formulae [J, and 
its exponential decay in both methods is presented in Fig. 
4. In the case of the puncture method, we have exhibited 
the results for distinct, but fixed values of the radial trun¬ 
cation order at the first region, namely = 20,30. In 
the later case the convergence is better. 

For the sake of completeness, we have exhibited in In 
Fig. 5 the behavior of the ADM mass in function of 
the amplitude Aq for n = 4, r7o = lifl = l and c = 
I, — I, —2. We noticed the same counterintuitive behavior 
of the ADM mass found in Refs. ii, that is, Madm 
initially decreases when Aq increases for Aq > 0. 

The last numerical test is to show the convergence of 
the L2-norm of Res(x, y, (/>), i. e., the residual equation 
associated to the Hamiltonian constraint We have 
considered the inversion method and introduced the vari¬ 
ables X, y as before. We show the results in Fig. 6 
corresponding to two cases. In the first, we have fixed 
Ny = 10 and increase for some values of the ampli¬ 
tude Ao, namely Aq = 0.01,0.5,1.O. Notice that L2(Res) 
achieves a limit minimum value after some value of 
that depends on Aq. In the second plot we have explored 
this aspect by fixing Aq = 1 , and for each value of Ny, 
evaluate the norm choosing N^ such as to get the mini¬ 
mum value, according to the upper graph. For instance, 
as we can see from the upper graph, for Ny = 10, we 
set Nx = 13 (for Aq = I). In both cases we have set 
<J = r]o = c = I, n = 4, and the decay of the L2(Res) is 
indeed exponential. 

We now turn to the problem of gravitational waveform 
extraction, whose accurate calculation is the most rele¬ 
vant problem in numerical relativity. The spin-weighted 
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FIG. 4. First: L2((I4'), inversion method. We have fixed 
Ny = 6 and the difference corresponds to the approximate 
solutions with + 5 and N^', the error L 2 is evaluated for 
r > a. Second: L 2 {d'i>), = 20 (box), = 30 (cir¬ 

cle); puncture method with domain decomposition; Ny = 6 
remains fixed and the difference corresponds to the approxi¬ 
mate solutions with + 5 and N ^.; the error L 2 is evaluated 
for r > a, i.e., the second domain. In all cases Lq = 9.0. 





FIG. 5. ADM masses evaluated for the first initial data with 
c = —2, —1,1, respectively, from top to bottom. 


scalar 'I'4 defined in the Newman-Penrose formalism [s^ 
provides a measure of the outgoing gravitational radi¬ 
ation [33, [3^. Therefore, we shall determine the pat¬ 
tern of radiation perceived by a distant observer from 
the source by examining the dominant terms resulting 
the limit r —^ 00 of 'I'4. This scalar is expressed by the 
following projection of the Weyl tensor C^^ap- 

= (30) 

where and m'' belong to the null tetrad basis adopted 
by Bernstein et al 0 and shown in the Appendix. The 
complete expressions for the real and imaginary parts 
of 'I'4 in the initial slice are also shown in the Appendix. 
Taking into account the asymptotic expression of the con¬ 
formal factor, the choices of unit lapse and zero shift, and 
since the function q(r, 9, /p) decays exponentially with r, 
we have found that d>4 ^ 0{r~^) at the initial slice, in¬ 
stead the typical decay 0{r~^) which characterizes the 
wave zone. We attribute this behavior to the condition 


of time symmetry demanding that Kjk = 0 at the initial 
slice. We may recover the standard asymptotic behav¬ 
ior of 'I'4 at subsequent slices in a dynamical setting. 
Notwithstanding this fact, we can consider the dominant 
terms of 1114 for large r as characterizing the structure of 
gravitational wave distribution at the initial slice. The 
asymptotic expression of 'I'4 is, 

lim r^\E'4 = -{r'i>),ee + (r'E') « cot(0) -I- (r'E') aa csc^ 9 + 

r—^oo 

3 ? 

+ ir'if),4,cot9]. (31) 

We have considered the approximated conformal fac¬ 
tor with Nx = 20, Ny = 16 corresponding to the first 
initial data and obtained the expressions for the above 
asymptotic real and imaginary parts of ^^4. Hereafter, 
these pieces are denoted by and dl™, respectively. 
We also have set n = 4, 770 = 1 and cr = 1 reducing 
the parameter space to (Ao,c). In general, the angular 
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FIG. 6. Exponential decay of the L 2 -norm associated to the 
residual equation in two situations. In the first, we have fixed 
Ny = 10 and increased for three values of the amplitude, 
^0 = 1.0, 0.5, 0.01, that corresponds to the curves from up to 
down. The exponential decay of the L 2 -norm is also observed 
for fixed = 1 by increasing Ny and Nx until achieving the 
saturation value for the combination Nx, Ny. 


distribution of IE'4 has the same symmetric of the con¬ 
formal factor, but it depends on the parameters {Aq,c). 
For the sake of convenience, all graphs correspond to 
Aq = —1, and the parameter c can assume one of these 
values: —2, —1, 0,1. In Fig. 7 we show the three and bi- 
dimensional polar plots of dt4 for the axisymmetric case 
(c = 0). In Fig. 8 we present a sequence of three and 
bi-dimensional polar plots of It can be seen the 

role of the parameter c in changing the angular or lobe 
structure of the pattern. On the other hand, the angu¬ 
lar pattern of ^'4™, shown in Fig. 9, does not change 
significantly with respect to c. 


V. FINAL REMARKS 

In this work, we have implemented a Galerkin- 
Collocation spectral algorithm to solve the Hamiltonian 
constraint corresponding to three-dimensional distorted 
black holes. These configurations can describe two plan- 



3 7t 


FIG. 7. From top to bottom we show the three and bi- 
dimensional polar plots (plane <() = 0) of 4*4 for the axisym¬ 
metric case (c = 0). The lobe structure is symmetric as ex¬ 
pected from the time symmetry condition imposed on the 
initial data. Also, it can be seen that the quadrupole mode is 
dominant. 

sible astrophysical situations: the late stages of black 
hole coalescence or the interaction of a black hole with a 
cloud of gravitational waves. 

We have solved the Hamiltonian constraint in the 
realm of the inversion method as a direct generaliza¬ 
tion of the axisymmetric case Q, and also implemented 
the Galerkin-Collocation version of the puncture method. 
According to this method, the black hole interior is not 
excised, and the spatial domain 0 < r < 00 is taken 
into consideration. We have developed a simple version 
of the domain decomposition when com par ed with other 
versions found in the literature (THIOI fl3l - [lq . In this 
case, we have divided the spatial domain into two re¬ 
gions, 0 < r < ro and r > rg- The numerical experiments 
indicated that the boundary r = rg is better placed at 
the black hole throat. For the sake of comparison, we 
have exhibited the convergence of the ADM mass for 
both, puncture with domain decomposition and inver¬ 
sion methods, in Fig. 3. 

In both methods the conformal factor is expressed in 
as a series expansion of radial and angular basis func¬ 
tions constituted by, respectively, a suitable combination 
of the rational Ghebyshev polynomials and spherical har¬ 
monics. With respect to the angular basis, we, metaphor- 
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FIG. 8. Sequence of three and bi-dimensional polar plots (plane </> = 7r/2) with c = —2,—1,1, from left to right. In 

all cases Ao — — l,n = 4 ,770 = l,o- = 1 in the first initial data. The deviation from axisymmetry produces a rich multipole 
structure. 



FIG. 9. The structure of 4*™ does not change considerably 
with c, Aq. 


ically speaking, are dealing with a double-edged sword. 
These functions are the best basis on the sphere offering 
exponential convergence for regular functions defined on 
the sphere. On the other hand, spherical harmonics are 
more complicated than an y ot her basis functions because 
they are two-dimensional [27|. 


It is worth commenting the similarities and differences 
of the present numerical implementation with the one by 
Pfeiffer et al. Q since it is a spectral code with simi¬ 
lar basis functions. The first difference is that we have 
adopted radial basis functions that are suitable linear 
combinations of the pure Chebyshev polynomials, instead 
of pure Chebyshev polynomials of Ref. Q- The linear 
combination is such that each basis function satisfies de 
boundary conditions. As pointed out by Heinrichs [s^, 
a combination of pure Chebyshev polynomials produces 
accurate results due to lower accumulated round-off er¬ 
ror when solving higher order differential equations. An¬ 
other difference comes from the use of mappings. We 
have considered the algebraic map [13 that is more ad¬ 
equate to describe functions with algebraic asymptotic 
behavior, which is the case of the conformal factor. In 
spite of considering the same angular basis of Ref. Q, 
we have employed the Galerkin method with numerical 
integration, GN-I, in the angular domain instead of the 
Collocation method. However, the Collocation method 
was used in the radial domain. Finally, we point out 
that in Ref. Q the spatial domain is divided into several 
regions, whereas we have considered only two regions. 
The simplicity is due to the specific initial data problem 
we are dealing. 
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We remark that the conformal factor is defined in 
the entire spatial domain under consideration in each 
method. This feature provides a natural and simple de¬ 
termination of the ADM mass by calculating the asymp¬ 
totic spatial limit of the term /dr (cf. Eq. 28). In 
addition, we have examined the angular pattern associ¬ 
ated to the dominant term of the spin-weighted scalar 
^'4 at the spatial infinity. In the present case we have 
found that ^'4 ~ 0{r~^). We can interpret such pat¬ 
terns as the indicators of the gravitational waves at a 
large distance from the distorted black hole. The an¬ 
gular patterns present a rich structure that depend upon 
the parameters of the initial data, and can be understood 
as gravitational-wave fingerprints distorted black holes. 

The Galerkin-Collocation method is a viable alterna¬ 
tive to solve the initial data problem of distorted three- 
dimensional black holes. The next natural direction of 
the present research is to study the dynamics of dis¬ 
torted black holes. There are two issues we will be focus¬ 
ing, namely, the gravitational wave templates produced 
in this dynamical process and the efficiency of the grav¬ 
itational wave extraction. The previous works on the 
dynamics of distorted non-rotating [s^, and rotating 
black holes have not discussed in details these issues. 


Appendix 


We present here null tetrad basis assuming unit 
and zero shift on the initial slice, and the with the 
metric given by Eq. ([5]). Then, 
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